%% Code description
% This code estimates the cost parameters.
% Table 1, Table 2, Table A.4, and Table A.5 are generated. 

%% Housekeeping
clear;

% root directory: change accordingly
root_dir = '...';

% add the demand estimation folder
addpath(fullfile(root_dir, 'matlab/demand_estimation_AK'));

%% Load data & demand estimation outputs
temp = load(fullfile(root_dir, 'matlab/demand_estimation_AK/demand_estimates.mat'));

K = temp.K;
D = temp.D;
N = temp.N;

% predicted demand behaviors at the estimates
Delta_xi       = temp.Delta_xi;       
s1_predicted   = temp.s1_predicted;   
s1_predicted_y = temp.s1_predicted_y; 
s1_predicted_y_c = temp.s1_predicted_y_c; 
s1_predicted_y_m = temp.s1_predicted_y_m; 
pr_stay        = temp.pr_stay;       

% demand estimates
alpha     = temp.alpha;
lapse_u   = temp.lapse_u; 
variety_u = temp.variety_u;

% other parameters
myopic = temp.myopic;

% display
disp(' ')
disp('***************************************************')
disp('Supply estmation starts')
disp('***************************************************')
disp(['alpha = ', num2str(alpha)])
disp(['variety_u (now fixed) = ', num2str(variety_u)])
disp(['N.ist = ', num2str(N.ist)])
disp(['Mean n_insurers by rep. fringe  = ', num2str(mean(D.n_insurers(D.major==0)))])

clear temp;

%% Parameters that are fixed
max_p1 = max(D.p1);
max_p2 = max(D.eta_mm)*max(max(D.p2_actual)); 

[pr_fc, x_fc, lapse, no_insurer_fe,...
    beta_c, beta_f, T1, T2, B1_c, B2_c, B1_f, B2_f, ...
    crra, rho, sigma_c1, sigma_c2, eulerc, ...
    N_y, pr_y, y_medi, delta, resource_y1, resource_y2] ...
    = set_parameters(K, max_p1, max_p2);

% type distribution among the insured in the HRS from HRS_stats.do (tab
% type if rhiltc==1)
HRS_dist = [0.0731;0.1603;0.3516;0.0464;0.1077;0.2610 ];
HRS_dist = HRS_dist/sum(HRS_dist);
pr_fc_n = pr_fc./pr_fc(3);
risk_temp = sum(HRS_dist.*pr_fc_n);
cost_weight = (pr_fc_n./risk_temp);

% decide whether to have type-dependent claims
use_cost_weight = 1;

%% Data dimension change I: Expand aggregate states
% - for data inputs pertaining to the 2nd stage, columns change from K to K*M
% get data
p1 = D.p1; % N.ist x 1
p2 = D.p2_actual; % N.ist x K.
r2 = D.r2; % N.ist x K, rate increases in %
eta_mm = D.eta_mm; % 1 x M rate increase multipliers
gamma_mm = D.gamma_mm; % 1 x M claims multipliers
D_claims2 = D.claims2; % N.ist x K
D_pr_s = D.pr_s; % N.ist x K, pr dist of k 
D_pr_mm = D.pr_mm; % N.ist x M, pr dist of mm =(interest rate, NH cost)
M = length(eta_mm);
N_ist = N.ist;
N_st = N.st;

% expand cols to K*M
[r2_mm, p2_mm, delta_p_mm, s2_mkt_mm, claims2_mm, pr_s_mm] ...
    = shell_expand_agg(r2, p1, s1_predicted, D_claims2, D_pr_s, D_pr_mm, ...
    eta_mm, gamma_mm, N_ist, K, M, delta);

%% Data dimension change II: Unmerge fringes
% - rows change from N.ist to NN:= (N.ist-N_st) + sum(D_nF)
D_unmerge = shell_unmerge_fringe(D, N_st, ...
    Delta_xi, s1_predicted, s1_predicted_y, s1_predicted_y_c, s1_predicted_y_m,...
    r2_mm, p2_mm, delta_p_mm, s2_mkt_mm, claims2_mm, pr_s_mm);

%% Step 1: Estimate rate stability regulation cost by MLE
s2_mkt_input  = D_unmerge.s2_mkt;
delta_p_input = D_unmerge.delta_p;
in_regression_input = D_unmerge.in_regression;
reg_after = (D_unmerge.reg_year+1 < D_unmerge.year); % 0 if reg_year NA
mkt_input = D_unmerge.mkt;

fc_value = N_st; % allow fixed cost to vary across markets

[rs_estimate, cost2_1, cost_rs, cost2_0_mat] ...
    = rstability_estimation(s2_mkt_input, delta_p_input, ...
    in_regression_input, reg_after, fc_value, N_st, mkt_input);

% variations in estimates by aggregate risks
analyze_rstability(cost_rs, cost2_1, M, K, D_unmerge);

%% Step 2: Estimate initial rate regulation cost using FOC
lr_target = 0.8;

Ep2 = sum(D_unmerge.pr_s.*D_unmerge.p2,2);
Er = mean(Ep2./D_unmerge.p1); 

% data input: all have been expanded (where appropriate) & unmerged
D_input = D_unmerge;

% cost1_e, cost_ml: NN x 1
[cost1_e, cost_ml] = minloss_estimation(lr_target, cost2_1, ...
    D_input, ...
    resource_y1, resource_y2, pr_y, N_y, ...
    K, M,...
    myopic, delta, B1_c, B2_c, B1_f, B2_f, ...
    beta_c, beta_f, T1, alpha, rho, crra, sigma_c1, Er, cost_weight, use_cost_weight);

%% Firm total profits
D_input = D_unmerge;

NN = size(cost2_1,1);
KM = size(cost2_1,2);
weight_claim_cost = zeros(NN, KM); %  claims*s2
pr_stay_fixed2 = (1-unique(delta))*ones(N_y,KM);
for j=1:NN
    s1_ind1_jj = D_input.s1_predicted_y(:,j); % N_y x 1
    s2_ind1_jj = repmat(s1_ind1_jj, 1, KM).* pr_stay_fixed2; % N_y x KM
    for k=1:KM
        D_claims2_m_jj = D_input.claims2(j,k); % scalar
        weight_claim_cost(j,k) = sum(cost_weight.*repmat(D_claims2_m_jj,N_y,1).*s2_ind1_jj(:,k).*pr_y, 'all'); %sum (N_y, KM) at firm level
    end
end

[profit, lr] = profit_ftn(D_input, cost_rs, cost_ml,...
    beta_f, T1, B1_f, B2_f, weight_claim_cost, use_cost_weight);

%% Step 3: Estimate market-specific entry cost distribution
D_input = D_unmerge;
N_st = N.st; % # of markets
N_f = 145; % potential entrants

% if 1, we estimate market-specific mu while fixing sigma. if 0, we estimate sigma intead. 
estimate_mu = 1; 

[mu_f_all, sigma_f_all,...
    mean_cost_entry, std_cost_entry, entry_exitflag_vec]...
    = entry_cost_estimation(D_input, profit, N_st, N_f, estimate_mu);

%% Save
save('supply_estimates.mat')

if fc_value==N_st
    % save to be used as initial conditions when re-estimating
    save('rs_estimate_ini_N_st.mat', 'rs_estimate');
    
    % export the estimates and run cost_estimates_regression.do in Stata.
    % this produces inputs needed to generate Table A.8.
    export_cost_estimates(cost2_0_mat, cost2_1, cost1_e, N_st, D_unmerge, D);
end

%% Export supply estimates to .tex file 
% Table 1, Table 2, Table A.4, Table A.5 are generated. 
MAIN_shell_export_tex









